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The linear noise approximation (LNA) offers a simple means by which one can study intrinsic noise 
in monostable biochemical networks. Using simple physical arguments, we have recently introduced 
the slow-scale LNA (ssLNA) which is a reduced version of the LNA under conditions of timescale 
separation. In this paper, we present the first rigorous derivation of the ssLNA using the projection 
operator technique and show that the ssLNA follows uniquely from the standard LNA under the 
same conditions of timescale separation as those required for the deterministic quasi-steady state 
approximation. We also show that the large molecule number limit of several common stochastic 
model reduction techniques under timescale separation conditions constitutes a special case of the 
ssLNA. 



I. INTRODUCTION 

In the study of complex systems, it is common to in- 
voke assumptions under which the dimension (and hence 
complexity) of the system is reduced; such a strategy of- 
ten leads to relatively simple theories capable of exact an- 
alytical predictions, which offer insights typically lost in 
numerical approaches. This approach is particularly use- 
ful in the study of biochemical pathway dynamics which 
typically involve the interaction of hundreds or thousands 
of different chemical species [1]. Deterministic models of 
such systems are frequently simplified by invoking the 
presence of well-separated timescales [1, 2], in partic- 
ular by means of the quasi-steady-state approximation 
(QSSA) [3, 4]. However it is well appreciated nowadays 
that a discrete stochastic approach in terms of chemical 
master equations (CMEs) is more faithful in describing 
dynamics inside living cells since the number of molecules 
of many species is in the range of few tens to few thou- 
sands [5]. The development of reduced stochastic mod- 
els consistently coarse-grained over timescales presents a 
significantly larger challenge than the reduction of de- 
terministic models because of the much larger number 
of differential equations which need to be solved in the 
former compared to the latter. 

Indeed it has been shown that it is not possible to 
formulate a reduced CME for the whole region of pa- 
rameter space in which the deterministic QSSA is valid 
but rather only over a subset of this space [6-8]. For 
example consider the Michaelis-Menten mechanism in 
which substrate reversibly binds with enzyme to form an 
enzyme-substrate complex which then decays into prod- 
uct. We are interested in conditions such that the sub- 
strate species decays over a much longer timescale than 
the enzyme and complex species. For this case, a consis- 



tently reduced CME can only be written provided that 
the complex decay rate into product is much less than its 
decay rate into substrate and enzyme which is not nec- 
essarily the case [8, 9]. However it is possible to write 
down reduced deterministic rate equations (REs) using 
the QSSA even when the aforementioned condition is not 
satisfied! [8] 

Recently we have shown that a reduced stochastic de- 
scription with the same range of validity as the determin- 
istic QSSA is possible. This reduction is achieved by first 
approximating the CME by linear Langevin equations, 
an approximation called the linear noise approximation 
(LNA) and then integrating out the fast fluctuations such 
that one obtains a reduced set of linear Langevin equa- 
tions for the fluctuations in the slowly varying species 
only. The latter is the slow-scale LNA (ssLNA) [10]. 
The advantage of the ssLNA over the reduced CME ap- 
proach is that the former is valid over the same range of 
parameter space as the QSSA, a claim which has been 
numerically verified for a number of biochemical circuits 
including cooperative and non-cooperative enzyme reac- 
tions and gene networks with or without negative feed- 
back regulation [10]. The ssLNA, as it currently stands, 
has been derived by means of simple and intuitively clear 
physical arguments but is lacking a formal and rigorous 
derivation. In this article we provide such a derivation 
and also show the relationship of the ssLNA with other 
common methods of stochastic model reduction based on 
timescale separation assumptions. 

The paper is divided as follows. In Sec. II, we con- 
sider the CME of a general system of elementary chemical 
reactions with two characteristic and clearly separated 
timescales and reformulate the conventional system size 
expansion of the CME under such conditions. In Sec. Ill 
we use the leading order terms of the expansion to show 
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FIG. 1. Schematic illustrating the derivation of deterministic 
and stochastic quasi-steady state approximations and rapid- 
equilibrium approximations for a reaction network with slow 
and fast species. The deterministic approach based on rate 
equations uses Tikhonov's theorem while the stochastic ap- 
proach utilizes the projection operator method applied to the 
system size expansion derived from a modified van Kampen 
ansatz. 



that the deterministic limit of the CME under timescale 
separation conditions is formally the same as the conven- 
tional reduced rate equations obtained from the QSSA. 
The next to leading order terms of the expansion pro- 
vide us with a Fokker-Planck equation which we reduce 
to a simpler Fokker-Planck equation for the slow fluctu- 
ations only, by means of the projection operator formal- 
ism. The reduced Fokker-Planck equation is obtained 
in closed-form for all monostable reaction networks and 
is the same as the ssLNA obtained in [10]. Finally in 
Sec. IV we show how the large molecule number limit 
of conventional stochastic reduction methods leads to a 
special case of the ssLNA and that hence these methods 
provide results in a more restrictive range of parameter 
space than the ssLNA does. The methodology of our 
approach is summarized in Fig. 1. 



II. THE SYSTEM SIZE EXPANSION UNDER 
TIMESCALE SEPARATION CONDITIONS 

We consider a system that comprises N chemical 
species (i = 1, . . . , N) confined in a compartment of 
volume ft and assume that the species can interact via 
j = 1, . . . , R elementary chemical reactions 



N 

£ r ijXi 

i=l 



(1) 



where and the stoichiometric coefficients [11] 

and kj is the macroscopic reaction rate. The constraint 
Si s ij — ^ for ah j implies that the reactions are uni- 
molecular or bimolecular and hence elementary. The to- 
tal number of species N is assumed to comprise N s slow 
and Nf = N — N s fast species, respectively. For conve- 
nience, we stick to the convention that X\ to Xjy a denote 
the slow species, while X^ s +i to Xpf are reserved for the 
fast species. 

In what follows, matrices are denoted by underlined 
quantities and all vectors are column vectors. Let Uj 
denote the number of molecules of species i, then the 
probability P(n, t) to find the system in a particular state 
n = (m, n N ) T is determined by the CME [12, 13] 



dt 



1 )f j (n,n)P(n,t) (2) 



-sn is the stoichiometric matrix 



where SVj = (S)y = 
and Ef is a step operator, the action of which on some 
function of the absolute number of molecules results in 
the same function but with rii replaced by rii 4- z, and 
fj (n, ft) are the entries of the microscopic rate function 
vector. The product £lfj(ri,Q)dt represents the propen- 
sity function, which has the meaning of the probability 
that reaction j takes place in a small time interval dt. 

The CME is typically analytically intractable and 
hence a systematic approximation method is needed. The 
system size expansion as developed by van Kampen is 
one such technique [12]. The heart of the method is the 
ansatz that the molecular concentration of the CME can 
be written as: 



n 
ft 



= ^+ft" 1 /2 



1, 



(3) 



where 4> is vector of concentrations as given by the cor- 
responding REs and the second term on the right hand 
side is a stochastic term describing fluctuations about the 
concentrations of the REs. 

We now impose timescale separation conditions, i.e., 
we assume that the transients in concentrations of a 
group of species decays much slower than those of the 
remaining group of species. The first group of species 
we label as the slow species while the second are the fast 
species. The characteristic timescales of these two are 
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t s and Tf, respectively. The vectors of molecular pop- 
ulations, of macroscopic concentrations and of fluctua- 
tions can be divided into subpopulations of slow and fast 
species: n = (ri s ,rif), <f> — (<f> s ,<fif) = (T s x s ,TfXf) and 

1/2-, _l/2 



ff= %,nf) = it* 

then be written as: 
1 n. 



< M ,Tf' tf). The ansatz Eq. (3) can 



= x 8 + (fir s 



T s Q 

i|=x / + (Or / )" 1 /V / 



(4) 



We now use this ansatz to derive an expansion of the 
CME in powers of the small parameter il -1 / 2 valid 
in the case of large volumes or, equivalently, for large 
copy numbers of molecules. For convenience the com- 
ponents of the various vectors will be denoted as fol- 
lows: x s = (xi,...,x N3 ) T , x f = (x Na+ i,...,x N ) T , e s = 
(ei, £n s ) T and tf = (£n s +i, £n) T ■ The probabil- 
ity distribution of molecular populations P(n, t) is hence 
forth changed into the distribution of slow and fast fluc- 
tuations II(e* s , tf,t). In what follows we consider how the 
time derivative, the step operator and the microscopic 
rate function vector which compose the CME Eq. (2) 
transform under the proposed ansatz Eq. (4). 



Transformation of the time derivative 



where we have partitioned the stoichiometric matrix into 
two parts: (S s )ij — (S)y for 1 < i < N s (the sto- 
ichiometric matrix for the slow species) and (Sf),y — 
(S)i+jv S) j for 1 < i < Nf (the stoichiometric matrix for 
the fast species). Note that g denotes some general func- 
tion of the molecule numbers. Applying the ansatz Eq. 
(4) it follows that the above equation can be written as 



JV 



IK 

i=l 



7(ei, ...,ejv) 



g{e x - (SlT.yV^SJu, ...,e Ns - ({It.)- 1 ^,)^, 

e^ + i-(nr / )- 1 /a(S / ) 1 j,...,e J v-(nr / )- 1 / 2 (S / ) JVlJ ). 

(8) 

The right hand side of the above equation can be Taylor 
expanded from which it follows that the step operator 
for the j th reaction has the following formal expansion in 
powers of the inverse square root of the system volume: 



i=l 

~ (_i)k n -*/2 

k=l 



~^VjSf 



(9) 



Using the change of variable theorem the time deriva- 
tive can be written as 



0_ 
dt 



P(n,t) 



dt 



dts 
dt 



' f 



dtj_ 
dt 



n(r, ,?/,*), (5) 



where V s = (d/dei, d/deN s ) T and V/ = 
(d/de Ng+ i, d/de N ) T . The time derivatives at con- 
stant n, denoted by can be computed from the ansatz 
Eq. (4), leading to 



d_ 
dt 



^ dS ^u(t s ,t f ,t). 



f dt 



(6) 



Transformation of the step operator 
By the definition of the step operator we have 

N 



l[E7 S *g( ni ,. 



,n N ) 



g(ni - (S s )ij, ...,Un s - ( S s )n 3 ,j, 
n N s +i — ( S f)i,j, ■■■,n N — (Sf)isr fl j), 



(7) 



Transformation of the microscopic rate function vector 

We consider four general types of elementary reactions 
depending on the order j of the reaction, for which the 
microscopic rate functions have been rigorously derived 
from microscopic considerations [14, 15]: (i) fj(n, 0) = 
kj, for a zeroth-order reaction by which a species is in- 
put into a compartment; (ii) fj(n, 0) = kjUufl^ 1 , for 
a first-order unimolecular reaction describing the decay 
of some species u; (iii) fj(n,Sl) — kjn u (n u — 1)0~ 2 for a 
second-order bimolecular reaction between two molecules 
of the same species u; (iv) fj(n,il) = kjn u n v Cl , for a 
second-order bimolecular reaction between two molecules 
of different species, u and v. 

With each of these cases one can also associate a 
macroscopic rate function, /, i.e., the rate of reaction 
as given by the corresponding rate equations. For the 
four cases discussed above these are: (i)fj(x S: Xf) = kj- 
(ii) fj(x s ,Xf) = k j( p u ; (iii) fj(x s ,x f ) = k^ (iv) 
fj(x s ,x f ) = kj(f> u 4> v , where <fo = t s xi if 1 < i < N s 
and fa = TfXi if N s + 1 < i < N. 

Given the microscopic and macroscopic rate functions 
for the four elementary reactions, one can write the for- 
mer in terms of the latter, leading to the general result 

fj(n,Q) = fj(S s ,S f ) + (nr^/^lf^Xf)^ 

+ (nrfr^^f^x^tf + oin- 1 ), (io) 
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where Vj s and V^, denote the vectors of derivatives with 
respect to the components of the vectors x s and Xf, re- 
spectively. 

This formula can be easily verified by considering each 
of the four elementary reactions discussed above. For 
example for case of a second-order j th reaction involving 
a slow species u and a fast species v, we have 

fj(n u ,n v ,Q) = k i-^-£ = kj(T s TfX u x v + (£It s )~ 1/2 t s x 

=f j (x u ,x v ) + (Or s )-V2 gte^) £u+ 
(Or/)- 1 ^ df ^> Xv) ev + 0(Q-i), (11) 

where we used the definitions of the microscopic and 
macroscopic rate functions given above and the ansatz 
Eq.(4). 



The transformed CME 

Substituting Eqs. (6), (9) and (10) in the CME Eq. 
(2) and rescaling time by the slow timescale r = t/r s , we 
obtain the transformed CME: 



d 



n (e s ,e/,r) 



( n Ts) V2 v t ^ _ g 8 / (lsif/) j n( e - s ,e),T) 

&Tf) 1/2 Vj (^L - 7 S //(a?., x f )) e>, r) 



fi 1 



( 7 £/ +7 1/2 Ant +£ s ) n(e s ,e/,r) + 0(Q" 



l/2> 



(12) 



where we have introduced the nondimensional ratio of 
slow and fast timescales 
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Note that the order fl° is defined by the operators 



-Vjj / e / + -VjD / V / 



(13) 

(14a) 
(14b) 

(14c) 



where J s = S.(V ar ./ r ) T , Jy, = S f (V Sf f r ) T and 
J s/ = S^V^/T, J /s = ^/(V^/^) 7, as well as 
D s = S s diag(/)S^, Bf = S/diag^Sj and D s/ = 
D/s = S,diag(/)SJ. 



III. DERIVATION OF THE SLOW-SCALE 
LINEAR NOISE APPROXIMATION 

A. Deterministic QSSA 

The leading order terms, C^fi 1 / 2 ), of Eq. (12) describe 
the dynamics of macroscopic concentrations and is given 
by the coupled set of REs: 



dx s 

1 dxf 
7 dr 



S s f(x s ,Xf) , 
Sff(x s ,x f ) . 



(15a) 
(15b) 



The presence of timescale separation is reflected by 
the large parameter 7 diminishing the time derivative in 
Eq. (15b). Such a set of equations present a special case 
in singular perturbation theory, where Eqs. (15a) and 
(15b) for the slow and fast variables are typically referred 
to as the degenerate and adjoined systems, respectively 
[11]. Tikhonov's first theorem [16, 17] states that a sim- 
plification of the above equations under timescale separa- 
tion conditions is possible whenever certain requirements 
are met: i) the solutions of both the degenerate and ad- 
joined systems, Eqs. (15) are unique and their right-hand 
sides are continuous functions; ii) the root Xf — h(x s ,r) 
is the stable solution of the adjoined system; iii) the ini- 
tial values x/(t = 0) are in the domain of influence of 
the solution as in ii). Whenever these prerequisites are 
met, the solution of the full system (15) for x s tends to 
the solution of the reduced system 



dx s 
~cfr~ 



= S s f(x s , h(x s )) . 



(16) 



in the limit of timescale separation, i.e., 7 1 — > 0. Note 
that Xf = h[x s ) is the solution of S ff(x s ,Xf) = 0. 

These requirements are typically fulfilled for the bio- 
chemical networks of interest. This is since the chemi- 
cal transformation scheme (1) is formulated for elemen- 
tary reactions, which are bimolecular or simpler, the right 
hand sides of Eqs. (15) are continuous polynomial func- 
tions of the second order at most. For monostable net- 
works, the rate equations admit a single steady state 
which is the same for the full and the reduced REs, i.e., 
Eq. (15) and (16). It is therefore clear that all the so- 
lutions will tend to this state with time, quicker for fast 
variables and slower for the slow ones. 



B. Adiabatic elimination of stochastic variables 
using the projection operator formalism 

In the previous subsection, we reviewed how imposing 
timescale separation on the deterministic level leads to 
reduced time evolution equations for the concentrations 
of the slow species. A related question which is the main 
issue of this article is: what is the reduced Fokker-Planck 
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equation describing the fluctuations about the concentra- 
tions predicted by the deterministic QSSA? 

The Fokker-Planck equation describing the fluctua- 
tions in the concentrations of both fast and slow species 
is given by the O(O ) terms in Eq. (12); this has the 
form 



— IL(e s , e>,r) = ( 7 £/ + 7 1/2 A„t + £ s ) n(e a , t f , r) . 

(17) 

From the definitions of the operators Eq. (14), it can 
be seen that £/ acts only on the fast variables, C s acts 
only on the slow variables and £; nt acts only on both 
slow and fast variables. Hence the three operators de- 
scribe processes evolving on fast, slow, and intermediate 
timescales, respectively. The dimensionless parameter 7 
here weights the degree of timescale separation in the sys- 
tem. The fact that 7 is the same as for the deterministic 
QSSA implies that the conditions for timescale separa- 
tion of the stochastic variables, in the limit of large fl, 
are exactly the same as those for the REs. 

For well-separated timescales, i.e., the case 7 ^ 1, we 
are typically interested in the probability distribution of 
slow variables, n(e s ,r) = J d€fU(e s , e/, r). Projection 
operator methods have been found useful in facilitating 
the adiabatic elimination of fast variables from stochastic 
descriptions [18, 19]. Here we use one such method to rig- 
orously obtain a reduced Fokker-Planck equation for the 
fluctuations in the concentrations of the slow variables. 

The main idea behind the projection operator method 
is that one specifies the quasi steady-state probability dis- 
tribution 7r(ej) of fast fluctuations, which is determined 
by the infinite 7-limit of Eq. (17): 

£ /7 r(e» = 0. (18) 

The reduction is then carried out defining the operator 

(Vn)(e s ,e f ,r) = 7r(e» J de f IL(e s , e), r) = 7r(e f )H(e s ,r), 

(19) 

projecting the probability distribution n(e* s ,e},r) onto 
the distribution of fast fluctuations evaluated at steady- 
state. Note that the above definition satisfies the relation 
V 2 = V and hence V is indeed a projector. 

In what follows we use the forms we have derived for 
the operators in the Fokker-Planck Eq. (17), i.e., those 
given by Eq. (14), to deduce three properties of the pro- 
jection operator. Given these properties we then show 
how the projection operator applied to Eq. (17) leads to 
a reduced Fokker-Planck equation in the slow variables 
only. 



Properties of the projection operator 

In this subsection we will show that the following prop- 
erties hold 

VC S = C S V , (20a) 
VCf=CfV = 0, (20b) 

vc int r = o. (20c) 

First, property (20a) follows from the fact that the 
projection operator V, as defined by Eq. (19), acts only 
on the fast variables e}, whereas C s , see Eq. (14a), acts 
only on the slow variables e s and hence the two operators 
V and C s commute. 

Second, we show that both equalities of property (20b) 
are satisfied. Considering the left hand side, we obtain 

VC f = n(e f ) J de)V / T (...) = 0, (21) 

since £/, as given by Eq. (14b), has the form of a di- 
vergence in the fast variables and hence by the partial 
integration lemma, its integral vanishes in the absence of 
boundary terms. By considering the right hand side, we 
have 

C f V=(C f 7r(e f ))Jde f =0, (22) 

by the quasi-steady state condition, Cfir = 0. 

The third property (20c) can be obtained as follows. 
The first, second and fourth terms of £- ln t as given by Eq. 
(14c) have the form of a divergence in the fast variables 
and hence by the partial integration lemma they give no 
contribution to 'P/CintP. The third term of Eq. (14c) also 
does not contribute albeit for a different reason than for 
the three other terms just discussed 

VC int V = n(e f ) J de } (-V T s l sS e f )-K{e } ) J de f 

= -n(e,Wl. f {?f)«J de f = 0. (23) 

Here, we have applied the partial integration lemma and 
then used the fact that {ef) n = 0. The latter follows 
from the explicit form of £/, see Eq. (14), which implies 
that the solution of the Fokker-Planck equation for the 
fast fluctuations Cfir(tf) = is a Gaussian distribution 
centered about zero. 



Derivation of the projection operator method 

We will now use the three properties of the projection 
operator just derived to obtain a reduced Fokker-Planck 
equation. Our approach in this subsection follows that 
of Gardiner [19, 20]. 



G 



We define Q = 1 — V and the following two functions 
v(t) = ra(e s ,e>,r), w(t) = QU(e s ,e f ,r), (24) 
together with their Laplace transforms 



v(s) 



dr e st v(t), w(s) 



dr e st w(t). 

(25) 



The latter has the distinct advantage that instead of deal- 
ing with differential equations we obtain a set of algebraic 
equations. Using Eqs. (24), (25) together with Eq. (17) 
we find 

sv(s) - v(0) = C s i{s) + -f 1/2 VC int w{s), (26a) 
sw(s) - w(0) = (C s +jC f + j 1/2 QC int )w(s) 

+ 7 1/2 QC int i(s). (26b) 

Note that use has been made of the properties (20a) and 
(20b). Solving for v(s), we obtain: 



sv(s) - v(0)-^ 2 VC int V( 7 )w(0) 

= [£ s + 7 ^Ant2%)£intHs), 



(27) 



where we have used definition (24) and property (20c) 
and introduced T>{^) = (s — C s — jCf — 7 1 / 2 <2£mt)~ 1 - 
From the above equation one can draw the limit 7—5-00 
for which T>(-y) ~ — (jCf) -1 : 



0_ 
dr 



v(t) 



v(t) 



(28) 



Using £; nt as given by Eq. (14) together with Eq. 
(29b) we can deduce the form of the reduced Fokker- 
Planck operator 



£ = L s -V T s l sf {e f C- f 1 V?)«l 



Vi l af (e f Cj L e/) v J s/ V s 



(30) 



Note that terms which have Vj to the left do not con- 
tribute to the reduced operator and hence are missing 
from the above equation. 

We proceed by evaluating the two distinct correlators 
appearing in the above expression explicitly. We shall 
make use of the identity 



due CfU = C 



l e C f u 



Cf(l-V), (31) 



which can be verified from straightforward integration 

and the fact that VH(el, e}, r) = lim^^ e CfU U(el, e), r) 

f T 
-f 



[19]. Using the fact that VeTix = 0, we can write 



(e) CJ 1 e f T )„ = J de f e f Cj\l - V)e f T n 
= ~ J du J de f e) e 

du(e f (u)e f T (0)) v 



Cru ?T _ 

e f 7T 



(32) 



where we have inverted the Laplace transform Eq. (25). 
Note that due to the vanishing of the third term on the 
left hand side of Eq. (27) this asymptotic limit is Marko- 
vian and hence does not require the knowledge of the 
initial distribution w(0) of the fast fluctuations. Using 
v(t) — 7r(e/)II(e^,r) and integrating over the fast fluctu- 
ations et we obtain 



an(e B ,T) 

dr 



C' = C s -(L mt CfC it 



(29a) 
(29b) 



where the angled brackets with subscript tt in Eq. (29b) 
denote the statistical average over the steady-state prob- 
ability distribution 7r(e/) of fast fluctuations. 



Derivation of the slow-scale linear noise approximation 



Note that in the third step we have taken into ac- 
count that ej T 7r(e/,u) = e Cj:u ef T Tr(ef,0) is a solution to 
d u ir = Cfir with the initial condition ejFn(0). One can 
utilize the Fourier transform of the autocorrelation ma- 
trix (e f (u)ef(0)) n = J duj/(2n)e lulu Sf(uj) to find that 



~S/(0) 



Similarly, one can show that 



1 _ _X _ T 



CJ 1 Vj), = - / due 3 -i u (e f V T f )« = — J 
Jo 



x -1 



(33) 



(34) 



Here, (e/Vj) w = — I with identity matrix I is evaluated 
by partial integration. Plugging Eqs. (33) and (34) into 
Eq. (30) gives a closed-form expression for the Fokker- 
Planck operator of the reduced probability distribution 
n(e„r): 



The above equation is a generic form for the Fokker- 
Planck equation for the slow fluctuations e* s under 
timescale separation conditions. What remains is to 
explicitly evaluate the average over 7r(e*y) such that we 
obtain a closed-form expression for the reduced Fokker- 
Planck equation. We now show these evaluation steps in 
detail. 



J =J 
D.. =D 



' ' ^ j T 



Usfl, 



% -1 ~ 



D /S f. 



(35a) 
(35b) 

(35c) 
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We complete our analysis by changing to natural vari- 
ables of concentration fluctuations ff s — (Ts/fi) 1 / 2 ^, of 
real time t = tt s and concentrations <fi s j = T s tx 8 j to 
finally yield the reduced Fokker-Planck equation for the 
fluctuations in the slow species 



J J s J s/ J j J /s ) 



(36a) 
(36b) 



- J S/ J7 1 D /S -(J S/ J7 1 D /S ) T , (36c) 

where V s denotes the derivative with respect to fj s . 
The coefficient matrices in the above expressions are 
given by J s = S,(V&/ r ) T . J / = S/^/T an d 
l Bf = S.<y$ f f r ) T , If. = Sf(y$ M f*) T as well as 
D s = Q^SsF Sj, D / = fi^S/F SJ, D s/ = Dj s = 

^ -1 S S FSJ and F = diag(/). 

We note that the slow-scale Jacobian (36b) coincides 
with the reduced Jacobian as obtained from the macro- 
scopic QSSA equations, i.e., Eqs. (f6) as shown in Ref. 
[10]. It is also important to note that the reduced diffu- 
sion matrix D ss admits the representation 



D. 



n-\A - b)(a - b; 



(37) 



where A = S s vT and B = J sf J j 1 S fy /Y. From 
this representation it can be immediately deduced that 
the reduced matrix D ss is symmetric and positive semi- 
definite, two crucial properties of the diffusion matri- 
ces for all Fokker-Planck equations [12]. Using standard 
methods [19] one can also obtain the Langevin equations 
corresponding to the slow-scale Fokker-Planck equation 
(36a). These are given by 

j t Ut) = Ji(i) + £T 1/2 (S S - J a/ J 7 X S /)>/!"?(*), 

(38) 

where the R dimensional vector T(t) is white Gaus- 
sian noise defined by (r^i)) = and {T i (t)r j :(t')) = 
Sij5(t — t'). Equations (36a) and (38) constitute equiva- 
lent forms of the LNA under timescale separation condi- 
tions, with the latter being particularly useful for Monte 
Carlo simulation purposes. This reduced LNA is pre- 
cisely the ssLNA introduced in Ref. [10]. 



which there exists slow and fast species whose transients 
decay on well-separated timescales, the rapid-equilibrium 
approximation divides the set of reactions into groups of 
slow and fast reactions where the latter determines the 
equilibrium of the fast species alone [11]. 

In the stochastic case there exist a variety of meth- 
ods separating fast and slow reactions. Examples are the 
"nested stochastic simulation algorithm" [23] , "slow-scale 
stochastic simulation algorithm" [24], the "stochastic 
partial equilibrium assumption" [25] and the "quasiequi- 
librium approximation" [26]. Confusingly, also the 
"stochastic quasi-steady state assumption" of Rao and 
Arkin [27] is valid only in the limit of slow and fast re- 
actions [8, 9]. All of these approaches have in common 
that the microscopic rate functions, or equivalently the 
propensities are rearranged into two groups: those asso- 
ciated with R s slow reactions which occur rarely over a 
long period of time and those associated with (R — R s ) 
fast reactions which occur frequently over a short period 
of time. Then we can define a constant /i such that 



max(/i,/ 2 ,...,/flJ 
min(/fl s +i,/iia+2, /r) 



» 1, 



(39) 



holds. For the sake of this article we will refer to approx- 
imations utilizing the above criterion as "rapid equilib- 
rium approximations" . Since our present derivation of 
the ssLNA is based only on the assumption of the pres- 
ence of slow and fast species it can be used to investigate 
the latter approximation as a partial case. 

Using the size parameter fi we can write the vector of 
macroscopic rate functions as 



(fs( 



x Sl x f ),nf f (x s ,Xf) 



(40) 



where f»(x s ,Xf) = (fi, fa, /rJ are the rates of the 

slow reactions and ff(£ a ,£f) = n~ 1 (f Rs+ i, f Rs+2 ,..., /«) 
are the rates of the fast reactions rescaled by the size 
parameter fi. Note that here fj, is determined by the 
infinite volume limit of Eq. (39). We can now partition 
the stoichiometric matrix into block matrices 



S = 



S/ 



qO) C(/) 
qW q(/) 



(41) 



IV. REACTION NETWORKS WITH SLOW 
AND FAST REACTIONS 

Besides the QSSA there exists another popular method 
to eliminate the fast variables from the determinis- 
tic REs; this is the rapid-equilibrium approximation 
[11, 21, 22]. While the QSSA considers the situation in 



discriminating slow and fast reactions (superscript) as 
well as slow and fast species (subscript). The matri- 
ces S i s ' and S denote the stoichiometrics of the slow 
and fast species in the slow and fast reactions, respec- 
tively, while Sg^ and S^ represent the stoichiometry 
of slow species in fast reactions and the stoichiometry of 
fast species in slow reactions. 



8 



A. Deterministic rapid-equilibrium approximation 

The macroscopic elimination starts from the conven- 
tional REs 

^ =MfiFV7+ £<*>/„ (42a) 

*f = ,s<«/} + s?>/„ (42b) 

which are similar to Eqs. (15) but discriminate slow and 
fast reactions by the size parameter \i. It is clear that 
fi and 7 must be of the same order since slow and fast 
timescales are determined by the size of the rate func- 
tions. 

Further, we observe that for // — > oo the above 
equations are not immediately of the form required by 
Tikhonov's theorem (compare Eqs. (15)), and hence the 
adiabatic elimination is not immediately applicable un- 
less we impose S i =0. This condition implies that the 
populations of slow species are not changed in fast reac- 
tions and is also imposed throughout the literature in re- 
ducing stochastic slow-fast reaction networks [24-26, 28]. 
Setting the time-derivative of the second equation to zero 
we can solve S / ff ~ for Xf = h(x s ) + Oip^ 1 ) to ob- 
tain the reduced system 

^=S^f s (x s ,h(x s )). (43) 

In the case where the equilibrium of the fast reactions 
S y // = is detailed balanced, the above approximation 
is called the deterministic rapid-equilibrium approxima- 
tion [11], i.e., the case when the fast reactions are given 
by a set of reversible reactions for which the forward and 
backward rates of each reaction cancel each other. 



It can be shown that the above expression coincides with 
the Jacobian obtained using Eq. (43). Third, we calcu- 
late the coefficients of the noise from Eq. (37) as 

A = S( s) VT7, B=0, (45) 

where F s = diag(/ s ) and Ff = diag(//). Note that 
the second equation follows from inserting S / y/TT — 

S { f s) y/Y~s + into th e definition of B after 

Eq. (37) together with the expression for Jj and J s f 
and taking the limit /i — > oo. This implies that B is of 
order /z _1//2 and hence the noise stemming from the fast 
reactions can be neglected in the limit /i —¥ oo. Finally, 
we can formulate the Langevin equations 

J t Ut)= + n- 1/a SW>/F;f(t). (46) 

Note that these are consistent with those obtained using 
stoichiometry S i and propensity vector f s (x s , h(x s )) of 
the reduced macroscopic Eqs. (43). Note that while in 
the ssLNA, Eq. (38), (which is consistent with QSSA 
conditions) both the noise in the fast and slow reactions 
contribute to the noise of the slow species, in the stochas- 
tic rapid-equilibrium approximation, Eq. (46), the noise 
in the slow reactions alone determines that in the slow 
species. 

The latter Langevin equation has also been obtained 
by Pahlajani et al. [28] starting from the decomposition 
of reactions into slow and fast categories. However our 
derivation is the first to clearly show that this Langevin 
equation is a partial case of the ssLNA and hence is only 
valid over a subset of the parameter space over which the 
QSSA holds. 

V. DISCUSSION 



B. Stochastic rapid-equilibrium approximation 

We can now apply the ssLNA to obtain the contri- 
bution of the fluctuations using the same conditions 
as used above for the deterministic rapid-equilibrium 
approximation. First, we make use of the condition 
Si f) = to obtain the coefficients 3 S = S { s s) (V $ t f?) T \ 
l f = H&f& h ff) T + SW(V ?/ /T) T and J s/ = 

S.P<y f( Z)\ If. = S«(V^/T) T + MSf (V $ JJF 
which distinguish contributions from slow and fast re- 
actions. Second, we can use these together within Eq. 
(36b) to obtain the reduced Jacobian by taking the limit 

fi — > oo 



= s W(v^/T) T 
sW(v^/T) T 



(44) 



In this article we have shown how to rigorously re- 
duce the linear noise approximation of the CME by using 
the projection operator formalism to eliminate the fast 
fluctuation variables. The resulting Langevin equation, 
Eq. (38) , is in agreement with the ssLNA as has been pre- 
viously derived from intuitive arguments only [10]. The 
present derivation provides a rigorous basis by deriving 
the LNA from the system size expansion using a modified 
van Kampen ansatz, Eq. (4) which is applicable under 
conditions of timescale separation. The resulting REs, 
Eq. (15) and the Fokker-Planck equation (17) obtained 
by this approach are of a particular form which allows 
direct application of Tikhonov's theorem and the projec- 
tion operator method, respectively. Hence by this proce- 
dure it is guaranteed that the mesoscopic elimination of 
the fast fluctuations is valid in exactly the same limit as 
the macroscopic elimination of the concentrations of the 
fast species by the deterministic QSSA. 

For reaction networks composed of slow and fast reac- 
tions, conditions considered by the majority of available 
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stochastic model reduction methods employing timescale 
separation, we have shown that the limit of large volumes 
or molecule numbers the CME can be approximated by 
a Langevin equation, Eq. (46), which is a special case 
of the ssLNA. The main advantage of using the ssLNA 
over the various aforementioned methods [23-28] is that 
the ssLNA is valid over a larger parameter range than the 
latter methods. The path integral approach developed in 
[29] also enjoys this property. However the ssLNA enjoys 
the further advantage that it is available in closed form 
for any monostable reaction network and hence can be 
readily constructed from knowledge of the stoichiomet- 
ric matrix and deterministic rate equations. The disad- 
vantage of the ssLNA is that for pathways composed of 
some second order reactions, the ssLNA (and the LNA on 
which it is based) is valid only for large enough molecule 
numbers [30, 31]. This limitation can be lifted by consid- 
eration of higher order terms in the system size expan- 
sion; such calculations present more formidable analyti- 
cal challenges than encountered in the derivation of the 
ssLNA and are under current investigation. 

For realistic biochemical networks there may be partic- 
ular parameter ranges for which the stability of the dy- 
namics is either monostable or bistable or even oscillatory 
states can occur. The applicability of the method there- 
fore depends on the type of stability realized for a par- 
ticular network under consideration. Oscillatory states 
are found for 10% of the transcriptome and 20% of the 
proteome in mouse liver [32, 33] and similar fractions in 
the human metabolome [34] . Although bimodal distribu- 
tions have been observed experimentally, as for instance 
in the lac operon of E. coli [35], a recent proteome-wide 
study suggests that such probability distributions (which 
potentially indicate bistability) are quite rare [36] and 



similarly for the human transcriptome [37]. Despite the 
fact that bistable and oscillatory properties are impor- 
tant for specific cellular functions, it appears that monos- 
table networks for which the present theory has been de- 
veloped are common in living cells. 

We note that the LNA has been applied also to net- 
works with limit cycles [38]. The resulting equation is 
still a linear Fokker-Planck equation and hence the elim- 
ination of fast variables can be performed along the same 
lines as in the present derivation. However, the anal- 
ysis of the resulting equation has to be carried out by 
means of Floquet theory due to the inherent phase dif- 
fusion in these systems [39]. For bistable systems, the 
underlying distribution cannot be captured by a linear 
Fokker-Planck Eq. and hence the LNA is not applicable 
in this case [12]. A commonly employed procedure to 
eliminate fast variables in such systems is the stochastic 
QSSA. A recent numerical study reported considerable 
discrepancies in the probability distributions of full and 
QSSA-reduced bistable systems [40]. We expect these 
discrepancies to be similar to the difference in the ssLNA 
and rapid-equilibrium approximations discussed in this 
article for the monostable case. Therefore, developing a 
technique to rigorously reduce the CME of bistable net- 
works remains still an open question. Since our method 
is devised for monostable systems and can be extended 
to oscillatory ones under timescale separation conditions 
we expect it to be of broad applicability for the study of 
intracellular reaction networks. 
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